library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(sdmTMB)
library(tidyr)

Mesh

# 300 knots
plot(mesh)

Model specifications

m1 <- sdmTMB(
  bio.adult ~ 0 + as.factor(Year) + GearCat + s(Month, bs = "cc", k = 10),
  data = dat,
  mesh = mesh,
  offset = log(dat$SweptArea),
  family = tweedie(link = "log"),
  spatial = "on",
  time = "Year",
  spatiotemporal = "IID"
)
## Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
## glmmTMB was built with TMB version 1.8.1
## Current TMB version is 1.9.0
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)
## Warning: The model may not have converged. Maximum final gradient:
## 0.0108123051405755.

Some diagnostics

dat$resids <- residuals(m1) # randomized quantile residuals

hist(dat$resids)

qqnorm(dat$resids[is.finite(dat$resids)])
qqline(dat$resids[is.finite(dat$resids)])

# check spatial autocorrelation
ggplot(filter(dat,Year%%5==1), aes(lon, lat, col = resids)) +
  scale_colour_gradient2() +
  geom_point() + geom_sf(data=map,inherit.aes = FALSE)+xlim(-12,30)+ylim(47,62)+
  facet_wrap(~Year,ncol=2) 
## Warning: Removed 758 rows containing missing values (geom_point).

Terms

m1
## Spatiotemporal model fit by ML ['sdmTMB']
## Formula: bio.adult ~ 0 + as.factor(Year) + GearCat + s(Month, bs = "cc", 
##  Formula:     k = 10)
## Mesh: mesh
## Time column: Year
## Data: dat
## Family: tweedie(link = 'log')
##  
##                     coef.est coef.se
## as.factor(Year)1970    11.34    1.24
## as.factor(Year)1971    10.15    0.69
## as.factor(Year)1972    10.10    0.67
## as.factor(Year)1973    10.01    0.83
## as.factor(Year)1974     9.99    1.43
## as.factor(Year)1975     9.08    0.82
## as.factor(Year)1976    10.04    0.74
## as.factor(Year)1977     9.53    0.71
## as.factor(Year)1978    11.58    0.52
## as.factor(Year)1979    11.61    0.51
## as.factor(Year)1980    12.24    0.45
## as.factor(Year)1981    13.35    0.47
## as.factor(Year)1982    11.87    0.46
## as.factor(Year)1983    12.14    0.45
## as.factor(Year)1984    12.45    0.44
## as.factor(Year)1985    11.03    0.43
## as.factor(Year)1986    11.23    0.42
## as.factor(Year)1987    10.88    0.42
## as.factor(Year)1988    11.42    0.41
## as.factor(Year)1989    12.21    0.41
## as.factor(Year)1990    10.85    0.40
## as.factor(Year)1991    10.92    0.40
## as.factor(Year)1992    10.67    0.40
## as.factor(Year)1993    11.15    0.39
## as.factor(Year)1994    10.88    0.38
## as.factor(Year)1995    10.50    0.38
## as.factor(Year)1996    10.38    0.38
## as.factor(Year)1997    10.63    0.37
## as.factor(Year)1998    10.63    0.38
## as.factor(Year)1999    10.60    0.37
## as.factor(Year)2000    10.32    0.37
## as.factor(Year)2001    10.08    0.36
## as.factor(Year)2002     9.61    0.37
## as.factor(Year)2003     9.65    0.35
## as.factor(Year)2004     9.71    0.36
## as.factor(Year)2005     9.46    0.36
## as.factor(Year)2006     9.57    0.36
## as.factor(Year)2007     9.93    0.36
## as.factor(Year)2008    10.16    0.35
## as.factor(Year)2009     9.47    0.36
## as.factor(Year)2010     9.92    0.35
## as.factor(Year)2011    10.17    0.34
## as.factor(Year)2012     9.84    0.35
## as.factor(Year)2013     9.72    0.35
## as.factor(Year)2014     9.80    0.35
## as.factor(Year)2015     9.81    0.35
## as.factor(Year)2016     9.93    0.35
## as.factor(Year)2017     9.59    0.35
## as.factor(Year)2018     9.16    0.36
## as.factor(Year)2019     8.59    0.36
## as.factor(Year)2020     8.47    0.37
## as.factor(Year)2021     8.32    0.36
## GearCatBT               0.25    0.18
## GearCatGOV_CL           1.54    0.17
## GearCatGOV_GG           0.97    0.22
## GearCatOTT             -3.26    0.53
## GearCatPEL             -2.05    0.75
## GearCatTV               0.89    0.23
## 
## Smooth terms:
##            Std. Dev.
## sds(Month)      0.43
## 
## Dispersion parameter: 665.08
## Tweedie p: 1.61
## Matern range: 1.35
## Spatial SD: 0.37
## Spatiotemporal SD: 0.26
## ML criterion at convergence: 248325.051
## 
## See ?tidy.sdmTMB to extract these values as a data frame.

Predict onto the grid

grid=grid %>% filter(!Year %in% c(1967,1968))
grid$Year=as.integer(grid$Year)
dat %>% group_by(GearCat) %>% summarise(n=n())
## # A tibble: 7 × 2
##   GearCat     n
##   <fct>   <int>
## 1 BAK_GG    520
## 2 BT      33421
## 3 GOV_CL  37009
## 4 GOV_GG   6353
## 5 OTT        87
## 6 PEL        88
## 7 TV      10896
grid$GearCat = 'GOV_CL'
# calculate the mode of month
# getmode=function(x){uniq=unique(x)
# uniq[which.max(tabulate(match(x,uniq)))]}
# getmode(dat$Month) Month set to 9
grid$Month = 9

p <- predict(m1, newdata = grid)

Spatial random fields

plot_map(p, "omega_s") +
  scale_fill_gradient2() +
  ggtitle("Spatial random effects only")
## Warning: Removed 4817 rows containing missing values (geom_raster).

Spatio-temporal random fields

plot_map(filter(p,Year%%5==1), "epsilon_st") +
  scale_fill_gradient2() +
  facet_wrap(~Year,ncol = 2) +
  ggtitle("Spatiotemporal random effects only")
## Warning: Removed 1169 rows containing missing values (geom_raster).

plot_map(filter(p,Year%%5==1), "exp(est)/1000") +
  scale_fill_viridis_c(
    trans = "pseudo_log",
    # trim extreme high values to make spatial variation more visible
    na.value = "yellow", limits = c(0, quantile(exp(p$est)/1000, 0.995)),'cod adult biomass \n (kg km^-2)'
  ) +
  ggtitle("Prediction (fixed effects + all random effects)",
    subtitle = paste("maximum estimated biomass density (whole time series) =", round(max(exp(p$est))/1000))
  )+facet_wrap(~Year,ncol=2)
## Warning: Removed 1169 rows containing missing values (geom_raster).